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Abstract 


A method for analysis of progressive failure in the Computational Structural 
Mechanics Testbed is presented in this report The relationship employed in this 
analysis describes the matrix crack damage and fiber fracture via kinematics-based 
volume-averaged damage variables. Damage accumulation during monotonic and 
cyclic loads is predicted by damage evolution laws for tensile load conditions. The 
implementation of this damage model required the development of two testbed 
processors . While this report concentrates on the theory and usage of these 

processors , a complete listing of all testbed processors and inputs that are required 
for this analysis are included. Sample calculations for laminates subjected to 
monotonic and cyclic loads were performed to illustrate the damage accumulation , 
stress redistribution , and changes to the global response that occurs during the 
loading history. Residual strength predictions made with this information compared 
favorably with experimental measurements. 


Introduction 

Laminated composite structures are susceptible 
to the development of microcracks during their 
operational lives. While these microcracks tend to 
aggregate in high stress regions and result in 
localized regions of reduced stiffness and strength, 
the microcracks can affect the global response of 
the structure. This change in the global structure in 
turn can create high stresses and increase damage 
accumulation in another part of the structure. Thus 
to accurately predict the structural response and 
residual strength of a laminated composite structure, 
the effects of the accumulating damage must be 
incorporated into the global analysis. The approach 
taken is to develop damage-dependent constitutive 
equations at the ply level. These equations are then 
employed in the development of the lamination 
equations from which the constitutive module of the 
structural analysis algorithm is constructed. This 
algorithm is executed in a stepwise manner in which 
the damage-dependent ply-level results are used in 
the calculation of the global response for the next 
load step. This report will describe two 
Computational Structural Mechanics (CSM) Testbed 
(COMET) processors that were developed for the 
performance of such an analysis. A brief review of 
the theory behind the processors is first presented. 
The usage of these processors is then demonstrated. 
Since this analysis requires the use of other COMET 
processors, this report serves as a supplement to the 
Computational Structural Mechanics Testbed User’s 
Manual (ref. 1). 

It should be noted that the current damage 
model capability, computer code version 1.1, is 
limited to matrix cracking and fiber fracture under 
tensile loads. 

Symbols and Abbreviations 

A laminate extensional stiffness matrix 

B laminate coupling stiffness matrix 


P mode I matrix cracking scale factor 

D laminate bending stiffness matrix 

dpara material parameter determined from 

experimental data 

E n lamina longitudinal modulus 

E 22 lamina transverse modulus 

F applied force 

G 12 lamina shear modulus 

y tensile fiber fracture scale factor 

k material parameter determined from 

experimental data 
N x applied load 

n material parameter determined from 

experimental data 

PDALC Progressive Damage Analysis of 

Laminated Composites 
R percent of maximum load 

On lamina longitudinal critical stress 

On lamina shear critical stress 

On lamina transverse critical stress 

\| / mode II matrix cracking scale factor 

tpiy ply thickness 

u longitudinal extension 

u°, v°, w° undamaged midplane displacements 

x, y, z displacement fields 

v 12 lamina Poisson’s ratio 
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Damage-Dependent Constitutive 
Relationship 

The damage-dependent constitutive relationship 
employed in the COMET analysis is based on a 
continuum damage mechanics model proposed by 
Allen, Harris, and Groves (refs. 2 and 3). Rather 
than explicitly modeling each matrix crack in the 
material, the averaged kinematic effects of the 
matrix cracks in a representative volume are 
modeled by internal state variables. These internal 
state variables for matrix cracks are defined by the 
volume-averaged dyadic product of the crack face 
displacement, u p and the crack face normal, n^ as 
proposed by Vakulenko and Kachanov (ref. 4): 



where a? u is the second order tensor of internal 
state variables, is the local representative 

volume in the deformed state, and S is the crack 
surface area. This product can be interpreted as 
additional strains incurred by the material as a 
result of the internal damage. From 

micromechanics it has been found that the effects 
of the matrix cracks can be introduced into the ply- 
level constitutive equation as follows (ref. 5): 

<7l = [<21 {££-«"}> (2) 

where are the locally averaged components of 
stress, [Q] is the stiffness matrix in ply coordinates, 

and {e^} are the locally averaged components of 
strain. The laminate constitutive relationships are 
obtained by integrating the ply constitutive 
equations through the thickness of the laminate to 
produce 

{AT} = [A] {eI} + [B]{KTz.} + {/"}, ( 3 ) 

{M} = [B]{el} + [D]{ KL } + {g M }, (4) 


material will produce midplane strain and curvature 
contributions equivalent to those resulting from the 
damage-induced compliance increase. 

As the matrix cracks accumulate in the composite, 
the corresponding internal state variables must 
evolve to reflect the new damage state. The rate of 
change of these internal state variables is governed 
by the damage evolutionary relationships. The 
damage state at any point in the loading history is 
thus determined by integrating the damage 
evolutionary laws. Based on the observation that the 
accumulation of matrix cracks during cyclic 
loading is related to the strain energy release rate G 
in a power law manner (ref. 8). Lo et al. (ref. 9) 
have proposed the following evolutionary 

relationship for the internal state variable 
corresponding to mode I (opening mode) matrix 
cracks: 

„ da? - - 

da^= — £*- kG n dN, (5) 

dS 

The term da? 22 reflects the changes in the internal 
state variable with respect to changes in the crack 
surfaces. This term can be calculated analytically 
from a relationship that describes the average crack 
surface displacements in the pure opening mode 

o 

(mode I) for a medium containing alternating 0 
o 

and 90 plies (ref. 5). The term G is the strain 
energy release rate calculated from the ply-level 
damage-dependent stresses. The material 

parameters, K and n , are phenomenological in 
nature and must be determined from experimental 

data (refs. 10 and 11). Because £ and n are 
assumed to be material parameters, the values 
determined from one laminate stacking sequence 
should be valid for other laminates as well. Since 
the interactions with the adjacent plies and damage 
sites are implicitly reflected in the calculation of the 
ply-level response through the laminate averaging 
process, equation 5 is not restricted to any particular 
laminate stacking sequence. 


where {N} and {M} are the resultant force and 
moment vectors, respectively; [A], [B], and [D] are 
the well known laminate extensional, coupling, and 
bending stiffness matrices, respectively (ref. 6); 
o 

{£l 1 is the midplane strain vector; and {k l } is 

the midplane curvature vector; {f^ } and {g M } 
are the damage resultant force and moment vectors 
for matrix cracking, respectively (ref. 7). The 

application of {f M ) and {g M } to the undamaged 


When the material is subjected to quasi-static 
(monotonic) loads, the incremental change of the 
internal state variable is assumed to be 


. A y . V) if £ki> EkUru ; 

[0 if £kl< Eklcm 


( 6 ) 


! I 


2 



where Eklcrit * s critical tensile failure strain and 

P,y, and \| / are scale factors that describe the load 
carrying capability of the material after the 
occurrence of mode I (opening mode) matrix 
cracking, fiber fracture, and mode II (shear mode) 
matrix cracking, respectively. The physical 
interpretation of equation (2) is as follows: as long 
as the strains in a material element (local volume 
element or finite element) are less than the critical 

strains, £ki C rit’ no damage exists and the internal 
state variables have a zero value. When the strains 
reach their critical value, the element is damaged 
and this damage is represented by an internal state 
variable whose value is proportional to the local 
strain. The proportionality is dependent on the 

scale factors f}, y, and \\f. Currently, the scale 
factors are chosen such that when fiber fracture, 
mode II matrix cracking, or mode I matrix cracking 
occur in a ply within an element, the longitudinal, 
shear, and transverse stresses for that ply in that 
element are 


a n=you 

(7) 

Ou = yOn 

(8) 

022 = ft 022, 

(9) 

where tfn , an, and CT 22 

are the lamina 


longitudinal, shear, and transverse critical stresses, 
respectively. 

Structural Analysis Formulation 

In order to simplify the formulation, it is expedient 
to consider the special case of symmetric laminates. 
With this assumption, the coupling stiffness matrix, 
[B], becomes the null matrix and the in-plane and 
out-of-plane laminate equations are decoupled. 
The laminate equations (3) and (4) are then 
substituted into the plate equilibrium equations to 
yield the following governing differential equations 
for the plate deformations: 


-Px = Ai 


d'u 

X ~dx> 

dV 


- n 2 „ 

a u 


d'u 


+ 2 A 16 - - + ^66 - 2 

axdy dy 

— o 

a v 


+ Al6 2 + (Al2 + A 66) 

OX 

df" df? 

+ Ax, z H 1 t" ■ 


dxdy 
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These governing differential equations are 
integrated against variations in the displacement 
components to produce a weak form of the damage 
dependent laminated plate equilibrium equations. 
By substituting the corresponding displacement 
interpolation functions into the weak form of the 
plate equilibrium equations, the following 
equilibrium equations in matrix form are produced 
(ref. 12 and 13) 


[K] {5} = { F A } + { F M } (13) 

where [K] is the element stiffness matrix, {5} is the 
displacement vector, { Fa } is the applied force 
vector, and { Fm ) is the damage-induced force 
vector resulting from matrix cracking. Note that 
the effects of the internal damage now appear on 
the right hand side of the equilibrium equations as 
damage-induced force vectors. 

Structural Analysis Scheme 

The progression of damage is predicted by an 
iterative and incremental procedure outlined in the 
flow chart shown in figure 1. The first block of 
figure 1 is a description of the finite element 
model. Block numbers are shown in the right hand 
corner of each box in the flowchart. Blocks 2 and 

3 are processors that calculate the element stiffness 
matrices and assemble and factor the global 
stiffness matrix. The compliance changes due to 
damage are accounted for by combining the 
damage induced force vector with the applied force 
vector and solving for the global displacements in 
equation (6). This solution process occurs in block 

4 and 5 and then the element stress resultants are 
computed in block 6. 
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18. For next load cycle, go to step 13; else stop. 


The ply-level strains and stresses are computed in 
block 7 and as long as the strains in a material 
element (local volume element or finite element) are 

less than the critical strains, £kl C rit’ no damage exists 
and the internal state variables have zero values. 
When the strains reach their critical value, the 
element is damaged and this damage is represented 
by an internal state variable. The damage state is 
updated and the ply-level stresses and strains are 
post-processed in block 8. The analysis iterates 
over blocks 4 through 8 until equilibrium is 
established, and then the next load step is applied. 

The implementation of this analysis into the 
COMET code can be accomplished with the 
development of processors DRF and DGI. These 
processors, as with other COMET processors, are 
semi-independent computational modules that 
perform a specific set of tasks. Processor DRF first 
calculates the damage resultant forces and moments 
and then incorporates them into the global force 
vectors. The second processor, DGI, post-processes 
the elemental stress resultants into ply level stresses 
and strains using the damage dependent constitutive 
relationship. With this information, the processor 
computes the damage evolution and updates the 
damage state for the next series of calculations. The 
remaining calculations can be performed with 
existing COMET processors. The following is a 
listing in order of COMET processor executions for 
this analysis: 

1 . Procedure ES defines element parameters. 

2. Processor TAB defines joint locations, 
constraints, reference frames. 

3. Processor AUS builds tables of material and 
section properties and applied forces. 

4. Processor LAU forms constitutive matrix. 

5. Processor ELD defines elements. 

6. Processor E initializes element datasets, creates 
element datasets. 

7. Procedure ES initializes element matrices. 

8. Procedure ES calculates element intrinsic 
stiffness matrices. 

9. Processor RSEQ resequences nodes for 

minimum total execution time. 

10. Processor TOPO forms maps to guide assembly 
and factorization of system matrices. 

1 1. Processor K assembles system stiffness matrix. 

12. Processor INV factors system stiffness matrix. 

13. Continue 

14. Processor DRF forms damage resultant force 
vectors. 

15. Processor SSOL solves for static displacements 

16. Procedure STRESS calculates element stress 
resultants. 

17. Processor DGI calculates ply level stresses and 
damage evolution. 


The usage and theory behind each of the 
existing processors can be found in The 
Computational Structural Mechanics Testbed User's 
Manual (ref. 1). The processors DRF and DGI are 
described in Appendices A and B of this report, 
respectively. With the exception of processor DRF 
and DGI, other processors from the COMET 
processor library can be substituted into the listing 
above to perform the tasks specified. 

Example Calculations 

Example calculations were conducted using 
COMET to illustrate the features of the progressive 
damage code. The first example demonstrates the 
effects of the evolving matrix damage on a crossply 
laminated composite plate subjected to constant 
amplitude fatigue loads. The dimensions and 
boundary conditions for the laminated plate are 
shown in figure 2. This plate was discretized into 
24 four-node quadrilateral EX47 shell elements 
(ref. 14). In this example, the plate has a [0/90] s 
laminate stacking sequence and the ply-level 
mechanical properties are listed in table L These 
properties corresponded to those measured for 
IM7/5260 (ref. 11). A maximum load of 2500 lb/in 
at an R-ratio of 0.1 was applied to the laminate. 
The COMET runstream and input, as well as a 
segment of the output, for this example can be 
found in the section entitled “Progressive Failure 
Analysis Input” of Appendix B. 

The predicted distribution of the mode I matrix 
M 

crack damage CL 22 * n the 90° plies is shown in 
figure 3. The damage was greatest at the narrow 
end of the plate since the component of stress 
normal to the fiber was highest in this region. The 
higher stresses further translated to a greater 
amount of energy available for the initiation and 
propagation of additional damage. This was 
reflected in the damage evolution along the length 
of the plate. However as damage accumulated in 

the plate, the stress gradient in the 90° plies became 
less steep as shown in figure 4. The similarity in 
stress resulted in relatively uniform changes to the 
damage state at the higher load cycles. For this 
laminate stacking sequence, the load shed by the 

damaged 90° plies was absorbed by the 0° plies. 
The consequence of this load redistribution is an 
increase in the global displacements as illustrated in 
figure 5. The redistribution of load to the adjacent 
plies will affect the interlaminar shear stresses as 
well. This could create favorable conditions for the 
propagation of delamination. 



The second example examines the effects of 
damage accumulation during cyclic fatigue loads 
on the stiffness of notched laminates. In this 
example, the notched laminates are tension fatigue 
loaded for 100,000 cycles. The notched (central 
circular hole) laminate is shown in figure 6. 
Symmetry was assumed about the length and width 
of the laminate so that only a quarter of the 
laminate was modeled by the finite element model . 
This model, also shown in figure 6, consisted of 
320 four-node quadrilateral EX47 shell elements. 
Two lamainate stacking sequences, a cross-ply 

[0/90d and a quasi-isotropic [0/±45/90] f , were 

considered. These laminates possessed the same 
ply-level material properties as the first example. 
(See table 1.) The maximum fatigue loads 
employed in sample calculations are listed in table 
2. The COMET runstream for the fatigue loading 
is similar to the previous example which is listed in 
Appendix B. In figure 7, the predicted stiffness 
loss for the open-hole geometry is compared to 
experimentally measured values of stiffness loss 
measured over a 4-in. gage length, symmetric about 
the open hole. 

The final example considers the effects of 
monotonically increasing loads on the residual 
strength of AS4/938, AS4/8553-40, and AS4/3501- 
6 notched laminates (ref. 15). The notched 
laminate shown in figure 8 is considered in this 
example. Symmetry was assumed about the length 
and width of the laminate so that only a quarter of 
the laminate was modeled by the finite element 
model. The four inch wide panels with one inch 
notches had a mesh (figure 8) that consisted of 278 
four node quadrilateral EX47 shell elements. All 
of the elements on the x=0 centerline from the 
notch-tip to the edge of the panel were the same 
size for the 4, 12, and 36 inch wide panels. The 
aminate stacking sequence was orthotropic 

+45/0/90/+ 30/0^. These laminates possessed 

the ply-level material properties shown in table 3 
(ref. 16 and 17). The applied load is incrementally 
increased with each load step to simulate a ramp up 
load input. Failure of the component is assumed to 
have occurred when the elements that span the 
width of the laminate have sustained a level of fiber 
fracture such that the analysis cannot reach 
equilibrium within a given load step. The load at 
which this condition exists is used to calculate the 
residual strength. The COMET runstream for this 
example is listed in Appendix C. The predicted 
residual strengths are shown along with 
experimental measurements in figure 9. 

The load redistribution due to damage 
progression is simulated using monotonic damage 
growth parameters (P, xp, and y). These parameters 


can be complicated algebraic functions describing 
the complex behavior of load redistribution due to 
matrix cracking and fiber fracture. However, for 
the purpose of developing the framework for this 
progressive damage methodology, simple constants 
were chosen for the parameters. So for this 
analysis, the damage growth law parameters govern 
the load redistribution in a way that is similar to the 
ply discount method. It is not reasonable to assume 
a 100% load redistribution at the instant of failure 
for mode I matrix cracking and tensile fiber 
fracture. Therefore, a 90% load redistribution was 
assumed, i.e. the local ply stress is only 10% of the 
critical ply strength (|J = 0.1 and y = 0.1). As the 
applied load increases, mode I matrix cracking and 
tensile fiber fracture internal state variables increase 
in proportion to the local strains. This results in a 
constant stress level (10% of the critical ply 
strength) in the damaged plies, illustrated in figure 
10, with the load redistributing to the surrounding 
plies and elements. 

Based on Iosipescu shear data (ref. 18) there is 
o 

a shear strain ( Yx2 ) where the behavior is no 
longer linear and becomes almost perfectly plastic. 
At this strain level, ip is equal to 1.0 to simulate 
elastic-perfectly plastic shear stress/strain behavior as 
illustrated in figure 10. This implies that as the 
applied load increases, the damaged ply carries 
100% of the critical shear strength while the 
additional stress transfers to the surrounding plies 
and elements. When the shear strain becomes 
cr 

catastrophic (y 12 )> ^ IS assumed to be equal to 0.1 . 
The load redistribution for shear is now similar to 
the case for mode I matrix cracking. 

Concluding Remarks 

This report describes a progressive failure 
analysis for laminated composites that can be 
performed using the Computational Structural 
Mechanics (CSM) Testbed (COMET) finite element 
code. The present analysis utilizes a constitutive 
model that describes the kinematics of the matrix 
cracks via volume averaged internal state variables. 
The evolution of these internal state variables is 
governed by an experimentally based damage 
evolutionary relationship. The nonlinearity of the 
constitutive relationship and of the damage 
accumulation process requires that this analysis be 
performed incrementally and iteratively. 

Two processors were developed to perform the 
necessary calculations associated with this 
constitutive model. In the analysis scheme, these 
processors were called upon to interact with existing 
COMET processors to perform the progressive 
failure analysis. This report, which serves as a guide 
for performing progressive failure analysis on 
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COMET, provides a brief background on the 
constitutive model and the analysis methodology in 
COMET. The description and usage of the two 
progressive failure processors can be found in the 
appendices of this report. 

The results from the example problems 
illustrated the stress redistribution that occurred 
during the accumulation of matrix cracks and fiber 
fracture. This in turn influenced the damage 
evolution characteristics, the global displacements, 


and the residual strengths. It should be noted that 
the current damage model capability is limited to 
mode I (opening mode) and mode II (shearing 
mode) matrix cracking and fiber fracture under 
tensile loading conditions. The inclusion of other 
damage modes such as delamination and 

compression failure mechanisms will provide a 
more complete picture of the failure process. 



Appendix A 

Processor DRF 

Al. General Description 

This processor calculates the damage resultant forces and moments caused by matrix cracking in 
laminated composites. These resultant forces and moments when applied to an undamage laminate will 
produce an equivalent amount of displacements and curvatures to those resulting from the matrix crack 
surface kinematics in a damaged laminate. This enables an analysis of the response of a damaged laminate 
without having to update the stiffness matrix each time the damage state changes. Matrix crack damage is 
modeled in this processor by volume averaged crack surface kinematics using internal state variables (refs. 2 
and 3). 

Processor DRF and processor DGI, which is described in appendix B, were developed to perform 
progressive failure analysis of quasi-static and fatigue loaded laminates in Computational Structural 
Mechanics (CSM) Testbed (COMET). Analyses from these processors are stored in two formats. One is in 
standard format that is accessed by opening the output file. The other is a data set, which is stored in a 
testbed data library, and provides data to processors and post-processors (ref. 1). In this analysis, processor 
DRF is used in conjunction with COMET analysis processors to determine the static displacement and 
elemental stress resultants for a laminated composite structure containing matrix crack damage. Processor 
DGI then calculates the damage-dependent ply stresses. The damage state is updated based on the ply 
stresses and this procedure is repeated for the next load cycle. 

Al.l. Damage Dependent Constitutive Relationship 

In this processor, the effects of the matrix cracks are introduced into the ply-level constitutive equations 
as follows (ref. 3): 


<7 l = [Q]{£l-CCl}, 


(Al) 


where CT^.j are the locally averaged components of stress, [Q] is the ply level reduced stiffness matrix, and 

M 

Ejjj are the locally averaged components of strain, a^j are the components of the strain-like internal state 

variable for matrix cracking and are defined 
by 


a ? ■= — f u. n. ds, 
L ,J y J i J 

V L 


(A2) 


where Vl is the volume of an arbitrarily chosen representative volume of ply thickness which is sufficiently 
M 

large that aLij do not depend on Vl, Uj are the crack opening displacements, nj are the components of the 
vector normal to the crack face, and S is the surface value of the volume Vl. The present form of the model 

assumes that an, an, and a 22 , the internal state variables representing the tensile fiber fracture, mode II 
(shear mode) matrix cracking, and mode I (opening mode) matrix cracking, respectively, are the only 
nonzero damage components. 


A1.2. Damage Dependent Laminate Equations 

The ply level strains are defined as follows: 



£u, = E'L-ZKl* 

(A3) 

£ Lyy £ Lyy ZK Lyy 

(A4) 

£ Lxy £ Lxy ^ Lxy 

(A5) 


where and Kl are the midplane strains and curvatures, respectively. The aforementioned ply strains are 
then substituted into equation (Al) to produce the ply-level stresses. Damage-dependent lamination 
equations are obtained by integrating these ply stresses through the thickness of the laminate (ref. 15). Next, 
the stiffness matrix in the laminate equation is inverted to produce 



A BViN-f**] 
B D\ 1 M-g M \ 


(A6) 


where [A], [B], and [D] are, respectively, the undamaged laminate extensional, coupling, and bending 
stiffness matrices. They are defined by the following equations from reference 6: 


r— *— i 

ju 

II 

■ K 
loT 

1 

C 1 

(A7) 

^ k~\ 

(A8) 

D k = \ 

(A9) 


where [Q ]k is the transformed reduced elastic modulus matrix for the k th ply in laminate coordinates. In 
equation (A6), N are the components of the resultant force per unit length and M are the components of the 

resultant moments per unit length. The variables {f^ 1 } and {g M ) represent the contribution to the resultant 
forces and moments from matrix cracking and are calculated from, 

{/"} = “X[e] (A10) 

*=1 

{*"} — i ){“'}, < A1I) 

L k = I 

where { 0t M } k contains the matrix cracking internal state variables for the k 1 * 1 ply. Thus given the forces N 
and moments M, as well as the damage variables in each ply, equation (A6) can be utilized to calculate the 
o 

midsurface strains £l and curvature Kl. 

A2. Processor SYNTAX 

This processor uses keywords and qualifiers along with the CLIP command syntax (ref. 1). Two 
keywords are recognized: SELECT and STOP. 


I IT 
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A2.1 Keyword SELECT 

This keyword uses the qualifiers listed below to control the processor execution. 


Qualifier 

Default 

Meaning 

LIBRARY 

i 

Input and output library. 

ELEMENT 

ALL 

Element type ( EX47, EX97) used in the analysis. Default is all 
element types found in LIBRARY. 

SREF 

1 

Stress reference frame. Stress resultants may have been computed 
in the element stress/strain reference frame (SREF=0) or in one of 
three alternate reference frames. For SREF=1, the stress/strain x- 
direction is coincident with the global y-direction. For SREF=3 the 
stress/strain jc-direction is coincident with the global ^-direction. 
Note that the processor currently must have the stress/strain 

coincident with the global jc-direction (SREF=T). 


PRINT 

i 

Print flag. May be 0, 1, or 2; 2 results in the most output. 

MEMORY 

2 000 000 

Maximum number of words to be allocated in blank common. 

This is an artificial cap on memory put in place so that the dynamic 
memory manager does not attempt to use all of the space available 
on the machine in use. 

DSTATUS 

1 

Damage state flag. If no damage, DSTATUS = 0. If matrix 
cracking (cyclic load), DSTATUS = 1. If matrix cracking 
(monotonic load), DSTATUS = 22222. 

XFACTOR 

0.0 

Increases the specified applied forces by this factor at every load 
step. This is used in the residual strength calculations. 


A2.2 Keyword STOP 

This keyword has no qualifiers. 

A3. Subprocessors and Commands 

Processor DRF does not have subprocessors 

A4. Processor Data Interface 

A4.1. Processor Input Datasets 

Several datasets, listed below, are used as input for processor DRF. 


Input Dataset 

Contents 

ELTS.NAME 

Element names 

OMB. DATA. 1.1 

Material properties including strain allowables 


9 

































LAM.OMB.*.* 

Laminate stacking sequence 

ES. SUMMARY 

Various element information 

PROP.BTAB.2.102 

ABD matrix 

WALL.PROP. 1 . 1 

Shell wall dataset 


Element directory dataset 

| DEF.xxxx.*.* 

Element definition (connectivity) dataset 


Internal state variable dataset 

xxxx.EFIL.*,* 

Element nodal coordinates and transformations 

APPL.FORC 

Applied force and moments at joints 

A4.2 Processor Output Datasets 

These Datasets used 

as output for processor DRF. 

Output Dataset 

Contents 

APPL.FORC 

Applied force and moments at joints 

mm 

Temporary damage resultant force dataset 


Damage resultant force dataset 


A5. Limitations 

Only EX47 and EX97 elements implemented with the generic element processor ESI will be processed 
by processor DRF. All other elements will be ignored. The stress reference frame must be coincident with 
the global ^-direction. 

A6. Error Messages 

Fatal errors will occur when any of the required datasets are missing from the input data library or when 
the stress resultants at the integration points are missing. 

Warning messages will be written and execution will continue when there is a missing or unreadable 
keyword or qualifier or if any of the original SPAR elements are encountered. 

A7. Usage Guidelines and Examples 

A7.1 Runstream Organization 

The following listing illustrates the organization of a progressive failure analysis that uses COMET. 
Because of the nonlinear nature of the damage-dependent constitutive equation, this analysis is performed in 
a stepwise manner. 

1 . Procedure ES Define element parameters. 

2. Processor TAB Define joint locations, constraints, reference frames. 

3. Processor AUS Build tables of material and section properties and applied forces. 

4. Processor LAU Form constitutive matrix. 




















5. Processor ELD Define elements. 

6. Processor E Initialize element datasets, create element datasets. 

7. Procedure ES Initialize element matrices. 

8. Procedure ES Calculate element intrinsic stiffness matrices, 

9. Processor RSEQ Resequence nodes for minimum total execution time. 

10. Processor TOPO Form maps to guide assembly and factorization of system matrices. 

1 1 . Processor K Assemble system stiffness matrix. 

12. Processor INV Factor system stiffness matrix. 

13. Continue 

14. Processor DRF Form damage resultant force vectors. 

15. Processor SSOL Solve for static displacements 

16. Procedure STRESS Calculate element stress resultants 

17. Processor DGI Calculate ply level stresses and damage evolution 

18. For next load cycle, go to step 13; else stop. 

A7.2. Progressive Failure Analysis Input and Output 

Please refer to Processor DGI for usage example (Appendix 1.2). 

A8. Structure of Datasets Unique to Processor DRF 
A8.1. DRFC.xxxx 

This data set is created by processor DRF and uses the SYSVEC format. See APPL.FORC.iset.l This 
dataset contains the damage resultant forces and moments corresponding to the given matrix cracking state. 

A8.2. DFCT.xxxx 

Data set DFCT.xxx is created by processor DRF and uses the SYSVEC format. See APPL.FORC.iset.l 
This dataset contains the damage resultant forces and moments from the previous load step and is used to 
restore the applied force vector to the initial value. 

A8.3. ISV.xxxx 

This data set contains the matrix cracking internal state variables at each layer. The xxxx is the element 
name. The data is stored in a record named ALP AM. 1. This record contains n items, where 

n = nlayer x nintgpt x nelt 

and nlayer is the number of layers in the model, nintgpt is the number of integration points for the element, 
and nelt is the number of elements. 

The data is stored in the following order: 

1 . a M j , internal state variable associated with fiber fracture. 

L< ii 

2. , internal state variable associated with mode I opening of the matrix crack. 

3. , internal state variable associated with mode II opening of the matrix crack. 

The data storage occurs for every layer, every integration point, and every element. 
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Appendix B 

Processor DGI 

Bl. General Description 

Processor DGI predicts the evolution of matrix crack damage in laminated composites for monotonical 
loads and cyclic fatigue loads. The processor also calculates fiber fracture under tensile load conditions. The 
matrix crack damage is represented in this processor by volume-averaged crack surface kinematics that use 
internal state variables (refs. 2 and 3). The evolution of these internal state variables is governed by a 
phenomenological growth law. 

This processor was designed to perform progressive failure analysis of laminated composite structures in 
the Computational Structural Mechanics (CSM) Testbed (COMET). At each load cycle, the elemental stress 
resultants for a laminated composite structure are obtained from COMET with the effects of matrix crack 
damage accounted for by Processor DRF. Processor DGI then postprocesses this information and uses the 
ply-level stresses to determine the evolution of matrix crack damage in each ply of the laminate. This 
procedure is repeated until the specified number of load cycles has been reached. 

BL1 Damage Dependent Constitutive Relationship 

In this processor, the effects of the matrix cracks are introduced into the ply-level constitutive equations 
as follows (ref. 5): 


M = [e]{e,-af} (Bl) 

where {ol} are the locally averaged components of stress, [Q] is the ply level reduced stiffness matrix, and 

{el} are the locally averaged components of strain. { 0Cl } are the components of the strain-like internal 
state variable for matrix cracking and are defined by 


n M _ 

W4 - 
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where Vl is the volume of an arbitrarily chosen representative volume of ply thickness that is sufficiently 
M 

large that a^j do not depend on Vl, u\ are the crack opening displacements, and nj are the components of 

the vector normal to the crack face. The present form of the model assumes that an,a? 2 >and # 22 , the 
internal state variables representing the tensile fiber fracture, mode II (shear mode) matrix cracking, and 
mode I (opening mode) matrix cracking, respectively, are the only nonzero damage components. 


For a uniaxially loaded medium containing alternating 0° and 90° plies, af 22 has been found from a 
micromechanics solution to be related to the far field normal force and crack spacing as follows (ref. 5): 
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where 
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p is the force per unit length that is applied normal to the fibers and 2t and 2a are the layer thickness and 
crack spacing, respectively. The C 2222 is the modulus in the direction transverse to the fibers and C 1212 is 
the in-plane shear modulus. Both moduli are the undamaged properties. 

B1.2. Damage Evolution Relationship 

Equation (B3) is used when the matrix crack spacing is known in each ply of the laminate. Since it is 
usually necessary to predict the damage accumulation and response for a given load history, damage 
evolutionary relationships must be utilized to determine the values of the internal state variables. The 

following relationship was used for the rate of change of the internal state variable a™ n in each ply during 
fatigue loading conditions (ref. 9): 
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M 

where d 0 CL 2 2 describes the change in the internal state variable for a given change in the crack surface areas, 

E and n are material parameters (refs. 10 land 11), N is the number of load cycles, and G is the damage- 
dependent strain energy release rate for the ply of interest and is calculated from the following equation: 

G = V t C w ,(e t „- <B 6 ) 

where Vl is the local volume. Interactions with the adjacent plies will result in ply strains £Lij» that are 
affected by the strains in adjacent plies. Thus, the strain energy release rate G in each ply will be implicitly 
reflected in the calculation of the ply-level response, so that equation (B5) is not restricted to a particular 
laminate stacking sequence. Substituting equation (B6) in equation (B5) and integrating the result in each 
ply over time gives the current damage state in each ply for any fatigue load history. 

When the material is subjected to monotonically increasing loads, the rate of change of the internal state 
M 

variable CC^ ls described by 
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for strains exceeding £ llm( , £ 2 2 ^ and7 12 , respectively. If none of the critical strains are exceeded, there is 

no damage. The (Xu is the updated internal state variable, af is the internal state variable for the previous 
damage state, and dan, daily and dan are the incremental changes in the internal state variables for tensile 
fiber fracture, mode I matrix cracking, and mode II matrix cracking, respectively. The monotonic damage 
growth parameters (P, \| /, and y) can be complicated algebraic functions describing the complex behavior of 
load redistribution due to matrix cracking and fiber fracture. However, for the purpose of developing the 
framework for this progressive damage methodology, simple constants were chosen for the parameters. So 
for this analysis, the damage growth law parameters govern the load redistribution in a way that is similar to 
the ply discount method. It is not reasonable to assume a 100% load redistribution at the instant of failure 
for mode I matrix cracking and tensile fiber fracture. Therefore, a 90% load redistribution was assumed, i.e. 
the local ply stress is only 10% of the critical ply strength (P = 0.1 and y = 0.1). As the applied load 
increases, mode I matrix cracking and tensile fiber fracture internal state variables increase in proportion to 
the local strains. This results in a constant stress level (10% of the critical ply strength) in the damaged plies, 
illustrated in figure 10, with the load redistributing to the surrounding plies and elements, 

o 

Based on Iosipescu shear data (ref. 16) there is a shear strain ( Y 12 ) where the behavior is no longer 
linear and becomes almost perfectly plastic. At this strain level, \| / is equal to 1.0 to simulate elastic-perfectly 
plastic shear stress/strain behavior as illustrated in figure 10. This implies that as the applied load increases, 
the damaged ply carries 100% of the critical shear strength while the additional stress transfers to the 

cr 

surrounding plies and elements. When the shear strain becomes catastrophic (Y 12 )» V iS assumed to be equal 
to 0.1. The load redistribution for shear is now similar to the case for mode I matrix cracking. 


B2. PROCESSOR SYNTAX 

This processor uses keywords and qualifiers along with the CLIP command syntax. Two keywords are 
recognized: SELECT and STOP. 

B2.1 Keyword SELECT 

This keyword uses the qualifiers listed below to control the processor execution. 


Qualifier 

Default 

Meaning 

LIBRARY 

1 

Input and output library. 

ELEMENT 

ALL 

Element type ( EX47, EX97) used in the analysis. Default is all 
element types found in LIBRARY. 

LOAD_SET 

1 

Load set; i of input data set STRS.xxxx.i.y. 

SREF 

1 

Stress reference frame. Stress resultants may have been computed 
in the element stress/strain reference frame (SREF=0) or in one of 
three alternate reference frames. For SREF=1, the stress/strain x- 
direction is coincident with the global y-direction. For SREF=3 the 
stress/strain ^-direction is coincident with the global z-direction. 
Note that the processor currently must have the stress/strain 
coincident with the global ^-direction (SREF=1). 

PRINT 

1 


STEP 

0 

Step number in nonlinear analysis (i.e., i in the STRS.xxxx.fO data 
set for nonlinear analysis). 



















MEMORY 

2 000 000 

Maximum number of words to be allocated in blank common. 

This is an artificial cap on memory put in place so that the dynamic 
memory manager does not attempt to use all of the space available 
on the machine in use. 

DSTATUS 

1 

Damage state flag. If no damage, DSTATUS = 0. If matrix 
cracking (cyclic load), DSTATUS = 1. If matrix cracking 
(monotonic load), DSTATUS = 22222. 

INCSIZE 

1.0 

Increment size used in damage growth law. 

NCYCLE 

1 

Cycle number 

NINCR 

1 

Increment number 


B2.2 Keyword STOP 

This keyword has no qualifiers. 


B3. Subprocessors and Commands 

None. Processor DGI does not have subprocessors. 


B4. Processor Data Interface 
B4.1 Processor Input Datasets 

Several datasets, listed below, are used as input for processor DGI. 


Input Dataset 

Contents 

ELTS.NAME 

Element names 

STRS .xxxx.i.j 

Element stress resultants. Record named INTEG_PTS must exist. 

OMB.DATA.l.l 

Material properties including strain allowables 

LAM.OMB.*.* 

Laminate stacking sequence 

ES .SUMMARY 

Various element information 

PROP.BT AB .2. 1 02 

ABD matrix 

WALL.PROP. 1 . 1 

Shell wall dataset 


Element directory dataset 


Element definition (connectivity) dataset 


Internal state variable dataset 

DGP.DATA.1.1 

Damage growth law parameters data set 


B4.2 Processor Output Datasets 


Output Dataset 


Contents 
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ISV.xxx.*.* 

Internal state variable data set 

PDAT.xxxx 

Ply-level stresses, strains, and damage state 


B5. Limitations 


Only EX47 and EX97 elements implemented with the generic element processor ESI will be processed 
by processor DGI. All other elements will be ignored. The stress reference frame must be coincident with 
the global ^-direction. 

B6. Error Messages 

Fatal errors will occur when any of the required datasets are missing from the input data library or when 
the stress resultants at the integration points are missing. (See section B4.1.) 

Warning messages will be written and execution will continue when there is a missing or unreadable 
keyword or qualifier or if any of the original SPAR elements are encountered. 

B7. Usage Guidelines and Examples 


B7.1 Organization of Progressive Damage Analysis on Testbed 


The organization of the COMET processors for a progressive failure analysis is shown below. The 
nonlinear nature of the damage-dependent constitutive equation requires that this analysis be performed in a 
stepwise manner. 


1 . Procedure ES 

2. Processor TAB 

3. Processor AUS 

4. Processor LAU 

5. Processor ELD 

6. Processor E 

7. Procedure ES 

8. Procedure ES 

9. Processor RSEQ 

10. Processor TOPO 

1 1 . Processor K 

1 2. Processor INV 

13. Continue 

14. Processor DRF 

15. Processor SSOL 

16. Procedure STRESS 

17. Processor DGI 

18. For next load cycle, 


Define element parameters. 

Define joint locations, constraints, reference frames. 

Build tables of material and section properties and applied forces. 
Form constitutive matrix. 

Define elements. 

Initialize element datasets, create element datasets. 

Initialize element matrices. 

Calculate element intrinsic stiffness matrices. 

Resequence nodes for minimum total execution time. 

Form maps to guide assembly and factorization of system matrices. 
Assemble system stiffness matrix. 

Factor system stiffness matrix. 

Form damage resultant force vectors. 

Solve for static displacements 
Calculate element stress resultants 
Calculate ply level stresses and damage evolution 
go to step 13; else stop. 


B7.2. Progressive Damage Analysis Input 

The following listing illustrates the input from a progressive failure analysis. The problem being solved 
is the uniaxially tensile loaded tapered laminated plate, figure 2, described in the main body of this report. 
The listing contains the main runstream plus a procedure file to perform the calculations for each load cycle. 


# 

# 

cp $CSM_PRC/proclib.gal proclib.gal .Copy procedure library 

chmod u+w proclib.gal 
testbed « \endinput 




♦set echo=off 

♦set plib=28 

♦open 28 proclib.ga! /old 

♦open/new 1 qoutput.101 


.Do not echo input 

.Open procedure library 
.Open output library 


tapered panel 

EX47 4 node quad elements 
24 nodes, 14 elements 


♦add pffc.clp .Add procedure for repeating calculations 

♦def/a es_name = EX47 .Element name 
*def/a es_proc = ES 1 .Element processor name 

♦call ES ( function= T DEFINE ELEMENTS’; es_proc = <es_proc> 
es_name=<es_name> ) 


fxqt TAB 

START 24 .24 nodes 

JOINT LOCATIONS .Enter joint locations 

1 0.0 0.0 0.0 20.0 2.0 0.0 8 1 3 

8 0.0 10.0 0.0 20.0 8.0 0.0 

CONSTRAINT DEFINITION 1 .Constraints: 

zero 1, 2,3,4, 5: 1,17,8 .Fixed end 

zero 6: 1,24 .Suppress drilling DOF 


[xqt AUS 

SYSVEC : appl fore 

1=1 : J=8 : 3750.0 
1=1 : J=16 : 7500.0 
1=1 : J=24 : 3750.0 


.Create input datasets 

.Applied forces 


TABLE(NI=16,NJ=1): OMB DATA 1 1 .Ply level material data 

IM7/5260 
1=1, 2,3, 4, 5 

J=l: 22.162E+6 0.333 1.262E+6 0.754E+6 0.754E+6 
I=6,7,8,9 

J=I: 0.754E+6 1.0E-4 1.0E-4 0.01 
1=10,11,12,13,14,15,16 
J=l: 0.0 0.0 0.0 0.014 0.014 0.014 0.0 

TABLE(NI=3,NJ=3,ITYPE=0): LAM OMB 1 1 .Section properties 

J=I: 1 0.006 0.0 
J=2: 1 0.012 90.0 
J=3: 1 0.006 0.0 

TABLE(NI=3,NJ= 1 JTYPE=0): DGPDATA 1 1 .Damage evolution data 

J=l: 1.1695 5.5109 3.8686E-7 


[xqt LAU 
ONLINE=2 


.Create constitutive matrix 


[xqt ELD .Define connectivity 

<es_expe_cmd> 

NSECT = 1 : SREF=1 


1 

2 

10 

9 

2 

3 

11 

10 

3 

4 

12 

11 

4 

5 

13 

12 

5 

6 

14 

13 

6 

7 

15 

14 

7 

8 

16 

15 

9 

10 

18 

17 

10 

11 

19 

18 

11 

12 

20 

19 

12 

13 

21 

20 

13 

14 

22 

21 

14 

15 

23 

22 
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15 


16 


24 


23 


[xqtE .Initialize element datasets 

stop 

♦call ES (function^ INITIALIZE') .initialize element matrices 

♦call ES (function='FORM STIFFNESS/MATL’) .Form stiffness matrices 

[xqt RSEQ .Resequence 

reset maxcon-12 

[xqtTOPO .Create maps 

[xqtK .Assemble global stiffness matrix 

rxqtINV .Factor the global stiffness matrix 

♦def/i ns_overwrite=<true> 

Call procedure to perform calculations at each cycle 

♦call PFFC ( es_proc=<es_proc> ; es_name=<es_name> ; — 

N_fcycl=l ; N_lcycl=2000 ; N_cylinc=5 ; - 
NPRT=100 ) 

♦pack 1 
[xqt exit 
\endinput 

B7.2.1. Procedure to perform the loop through the calculations for each load cycle 
(file name pffc.clp) 

♦procedure PFFC ( es„proc ; es_name ; — 

N_fcycl ; N_lcycl ; N_cylinc ; — 

NPRT ) 


N_fcycl: first fatigue cycle 
NJcycl: last fatigue cycle 
N_cylinc: cycle increment 
NPRT: output storage cycle increment 

begin loop here 

♦set echo=off 

*def icount = 0 .Initialize print counter 

♦DO :CYCLOOP SNCYL = <[N_fcycl]>, <[NJcycI]>, <[N_cylinc]> 
♦def icount = ( <icount> + 1 ) 

♦if < <icount> /eq <[NPRT]> /or <$NCYL> /eq 1 > /then 
♦def iprint = 1 
♦def icount = 0 
♦else 

♦def iprint = 0 
♦endif 

♦def delinc = <[N_cyIinc]> 

[xqt DRF .Calculate damage resultant forces 

select /PRINT = 0 
stop 

[xqtSSOL .Solve for static displacements 

Calculate elemental stress resultants 

♦call STRESS (direction=l; location^ INTEG_PTS; print=<false> ) 

♦if < <IPRINT> /eq 1 > /then 
[xqtVPRT .Print static displacements 

format = 4 
print ST AT DISP 
stop 


1 8 



[xqtDGI .Calculate ply level stresses, 

select /PRINT = 2 .strains, and damage evolution 

select /INC_SIZE = <delinc> 

select /N_CYCLE - <$NCYL> 

select /NINCR = 1 

select /NINCR = <$SNCYL> 

stop 

*endif 

♦if < <IPRINT> /ne 1 > /then 
[xqt DGI .Calculate ply level stresses, 

select /PRINT = 0 .strains, and damage evolution 

select /IN C_SIZE = <delinc> 

select /N_CYCLE = <$NCYL> 

select /NINCR = 1 

select /NINCR = <$SNCYL> 

stop 

♦endif 

♦if < <IPRINT> /eq 1 > /then .Store datasets for post processing 

♦copy 1, PLYDT.<ES_NAME>.<$NCYL>. 1 = 1, PDAT.<ES_NAME>.* 

♦copy 1, DISP.<ES_NAME>.<$NCYL>.1 = 1, STAT.DISP.* 

♦endif 

: CYCLOOP 

. *set echo=off 
♦end 


B7.3. Progressive Damage Analysis Output 

The following is a partial list of a progressive failure analysis output produced by processor DGI. Data 
for postprocessing is stored in dataset PLYDT.xxxx.xxx.l . 

♦♦BEGIN DGI ** DATA SPACE= 2000000 WORDS 
CYCLE NUM.= 496 


ELEMENT NUMBER 1 TYPE EX47 
EVALUATION (TNTG) POINT NUMBER 1 
REFERENCE SURFACE STRAINS AND CURVATURES 

E0-X E0-Y E0-XY K-X K-Y K-XY 

0.4619E-02 -0.6946E-04 0.1180E-02 O.OOOOE+OO 0.0000E+00 O.OOOOE+OO 
COMBINED BENDING AND MEMBRANE STRESSES, STRAINS, AND MATRIX CRACK DAMAGE 
VARIABLE 

FOR EACH LAYER OF ELEMENT 1 TYPE EX47 


LAYER 

THETA SIG-1 

SIG-2 

TAU-12 

STRAIN- 1 

STRAIN-2 GAMMA- 12 


1 

0 . 

0.103E+06 

0. 1 87E+04 

0.890E+03 

0.462E-02 

-0.695E-04 

0.118E-02 

2 

90. 

0.384E+03 

0.578E+04 

-0.890E+03 

-0.695E-04 

0.462E-02 

-0. 1 18E-02 

3 

0 . 

0.103E+06 

0. 1 87E+04 

0.890E+03 

0.462E-02 

-0.695E -04 

0.118E-02 


LAYER ALPM-11 ALPM-22 ALPM-12 

1 O.OOOE+OO 0.122E-1 1 0.000E+00 

2 O.OOOE+OO 0.473E-04 O.OOOE+OO 

3 0.000E+00 0.122E-1 1 O.OOOE+OO 


♦♦BEGIN DGI ** DATA SPACE= 2000000 WORDS 
CYCLE NUM.= 996 



ELEMENT NUMBER 1 TYPE EX47 

EVALUATION (TNTG) POINT NUMBER 1 

REFERENCE SURFACE STRAINS AND CURVATURES 

EO-X EOY EO-XY K-X K-Y K-XY 

0.4623E-02 -0.6882E-04 0.1183E-02 O.OOOOE+OO O.OOOOE+OO O.OOOOE+OO 

COMBINED BENDING AND MEMBRANE STRESSES, STRAINS, AND MATRIX CRACK DAMAGE 

VARIABLE 

FOR EACH LAYER OF ELEMENT 1 TYPE EX47 
LAYER THETA SIG-1 SIG-2 TAU-12 STRAIN- 1 STRAIN-2 GAMMA-12 

1 0. 0. 103E+06 0.187E+04 0.892E+03 0.462E-02 -0.688E-04 

2 90. 0.382E+03 0.573E+04 -0.892E+03 -0.688E-04 0.462E-02 

3 0. 0.103E+06 0.187E+04 0.892E+03 0.462E-02 -0.688E-04 

LAYER ALPM-II ALPM-22 ALPM-12 

1 O.OOOE+OO 0.246E-1 1 O.OOOE+OO 

2 O.OOOE+OO 0.901E-04 O.OOOE+OO 

3 O.OOOE+OO 0.246E-1 1 O.OOOE+OO 


** BEGIN DGI ** DATA SPACE= 2000000 WORDS 
CYCLE NUM.= 1496 


ELEMENT NUMBER 1 TYPE EX47 

EVALUATION (INTG) POINT NUMBER 1 

REFERENCE SURFACE STRAINS AND CURVATURES 

E0-X E0-Y EO-XY K-X K-Y K-XY 

0.4625E-02 -0.6839E-04 0.1184E-02 O.OOOOE+OO 0.0000E+00 0.0000E+00 

COMBINED BENDING AND MEMBRANE STRESSES, STRAINS, AND MATRIX CRACK DAMAGE 

VARIABLE 

FOR EACH LAYER OF ELEMENT 1 TYPE EX47 
LAYER THETA SIG-1 SIG-2 TAU-12 STRAIN- 1 STRAIN-2 GAMMA- 12 

1 0. 0.103E+06 0.187E+04 0.893E+03 0.463E-02 -0.684E-04 

2 90. 0.376E+03 0.568E+04 -0.893E+03 -0.684E-04 0.463E-02 

3 0. 0.103E+06 0.187E+04 0.893E+03 0.463E-02 -0.684E-04 

LAYER ALPM-11 ALPM-22 ALPM-12 

1 O.OOOE+OO 0.372E-1 1 O.OOOE+OO 

2 0.000E+00 0.129E-03 O.OOOE+OO 

3 O.OOOE+OO 0.372E-1 1 O.OOOE+OO 


** BEGIN DGI ** DATA SPACE- 2000000 WORDS 
CYCLE NUM.= 1996 


ELEMENT NUMBER 1 TYPE EX47 

EVALUATION (INTG) POINT NUMBER 1 

REFERENCE SURFACE STRAINS AND CURVATURES 

E0-X E0-Y EO-XY K-X K-Y K-XY 

0.4627E-02 -0.6806E-04 0.1185E-02 O.OOOOE+OO 0.0000E+00 O.OOOOE+OO 

COMBINED BENDING AND MEMBRANE STRESSES, STRAINS, AND MATRIX CRACK DAMAGE 

VARIABLE 

FOR EACH LAYER OF ELEMENT 1 TYPE EX47 


LAYER 

THETA SIG-1 SIG-2 TAU-12 

STRAIN- 1 

STRAIN-2 GAMMA- 12 


1 

0. 

0.103E+06 0.187E+04 

0.894E+03 

0.463E-02 

-0.681E-04 

0.119E-02 

2 

90. 

0.370E+03 0.564E+04 

-0.894E+03 

-0.681E-04 

0.463E-02 

-0. 1 19E-02 

3 

0. 

0.103E+06 0.187E+04 

0.894E+03 

0.463E-02 

-0.681E-04 

0. 1 19E-02 

LAYER 

ALPM-11 

ALPM-22 ALPM-12 





1 

O.OOOE+OO 

0.500E-1 1 0.000E+00 





2 

O.OOOE+OO 

0.164E-03 O.OOOE+OO 





3 

O.OOOE+OO 

0.500E-1 1 O.OOOE+OO 






T T I 


0.118E-02 
-0. 1 1 8E-02 
0.118E-02 


0. 1 18E-02 
-0.118E-02 
0.118E-02 


20 



B8. Structure of Datasets Unique to Processor DGI 


B8.1. PDAT .xxxx 


Data set PDAT. xxxx contains ply-level damage dependent stresses, strains, and internal state variables. 
Data are centroidal values. The variable xxxx is the element name. The data for each element is stored in a 
record named DAT_PLY.*Wr, where ielt is the element number. Each record contains n items, where 

n - nlayer x 9 

and nlayer is the number of layers in the model. 

The data is expressed with respect to ply coordinates and is stored in the following order: 

1. 0\\ normal stress in the fiber direction. 

2. O 22 normal stress transverse to the fibers. 

3. 0\2 shear stress. 

4. £] i strain in the fiber direction. 

5. £22 strain transverse to the fibers. 

6. £12 shearing strain. 

internal state variable associated with fiber fracture. 


M 

7 ■ a L1 1 
M 

a L22 
L 

«12 


8. 

9. 


internal state variable associated with mode I opening of the matrix crack, 
internal state variable associated with mode II opening of the matrix crack. 


Repeated nlayer times. 

B8.2. DGP. DATA. 1.1 


This data set is created by AUS/TABLE and contains the growth law parameters for the matrix cracking 
evolutionary relationship. The following variables are used to specify table size: 

NI = number of material parameters, for this case 3 
NJ = number of material systems, for this case 1 
Type = numerical format such as real or integer 

where NI and NJ are the number of columns and rows, respectively and Type specifies numerical format, real 
or integer. 

Each entry contains the following: 

1 . Growth law parameter, Tc . 

2. Growth law parameter, n . 

_ dotuj 

3. Parameter for determining , dpara 

These entries are repeated NJ times. 

B8.3. ISV.xxxx 


This data set contains the matrix cracking internal state variables at each layer. The variable xxxx is the 
element name. The data is stored in a record named ALPAM.l. 

This record contains n items, where 


n = nlayer x nintgpt x nelt 



and nlayer is the number of layers in the model, nintgpt is the number of integration points for the element, 
and nelt is the number of elements. 

The data is stored in the following order: 

M . 

1 . OCu ] , internal state variable associated with fiber fracture. 

M 

2. &122 > internal state variable associated with mode I opening of the matrix crack. 

M 

3. &L]2 > internal state variable associated with mode II opening of the matrix crack. 

The data storage occurs for every layer, every integration point, and every element. 



Appendix C 

Residual Strength Runstream 
Cl. General Description 

This appendix lists a sample runstream that was used to calculate the residual strength of 
an orthotropic center-notched laminate that was monotonically loaded to failure. 

C2. Residual Strength Analysis Input 

The following listing illustrates the input from a residual strength analysis. The problem being solved 
is the uniaxially tensile loaded center-notched orthotropic laminated plate shown in figure 8. The listing 
contains the main runstream plus a procedure file. The procedure file calculates the response during the 
monotonic loading to failure and is presented in this appendix. The finite element mesh was created using 
PATRAN. The file PT2T.PRC was created using the testbed PATRAN Tq Testbed (PT2T) neutral file 
converter. This file contains all of the nodal locations, connectivity matrix, boundary conditions, and 
applied forces. 


cp $CSM_PRC/proclib.gal proclib.gal 

chmod u-fw proclib.gal 

testbed > cct.o « \endinput 

♦set echo=off 

♦set plib=28 

♦open 28 proclib.gal /old 
♦open/new 1 cct.101 

rectangular panel with center-cut slit 
quarter panel mesh 
EX47 4 node quad elements 

*add pffdm.clp 
♦add initialize.clp 
♦ADD PT2T.PRC 
*def/a es_name = EX47 
*def/a es_proc = ES 1 

♦call ES ( function = DEFINE ELEMENTS'; es_proc = <es_proc> 
es_name=<es_name> ) 

[xqt TAB 
START 800 
♦call PT2TJLOC 

CONSTRAINT DEFINITION 1 
♦call PT2T_BC 

[xqtAUS 

SYSVEC : appl fore 
♦call PT2T_AF 

T ABLE(NI= 1 6,N J = 1 ) : OMB DATA 1 1 
AS4/938 
1 = 1 , 2 , 3 , 4, 5 

J=l: 19.60E+6 0.32 1.36E+6 0.72E+6 0.72E+6 
1=6,7, 8, 9 

J=1 : 0.72E+6 1.0E-4 1.0E-4 0.01 
1=10,11,12,13,14,15,16 
J=l: 0.0 0.0 0.0 0.0148 0.005 0.010 0.0 

TABLE(NI=3,NJ=13,itype=0): LAM OMB 1 1 
J=l: 1 0.0072 -45.0 
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J=2: 1 0.0072 45.0 
J=3: 1 0.0072 0.0 
J=4: 1 0.0072 90.0 
J=5: 1 0.0072 -30.0 
J=6: 1 0.0072 30.0 
J=7: 1 0.0072 0.0 
J=8: 1 0.0072 30.0 
J=9: 1 0.0072 -30.0 
J=10: 1 0.0072 90.0 
1=11:1 0.0072 0.0 
J=12: 1 0.0072 45.0 
J=13: 1 0.0072 -45.0 

[xqtLAU 
ONLINE=2 
[xqt ELD 

♦call PT2T_CONN 
♦def/a solver_name=INV 

♦call INITIALIZE ( renumber=l; rseq_method=-l; auto_dof_sup=0 ) 
♦call STIFFNESS ( type=MATL ) 

♦print 1 EX47.EFTL* DATA.l /m=100 

♦call FACTOR ( input_matrix=K; output_matrix=K ) 

*def/i ns_overwrite=<true> 

♦call PFFDM ( es_proc=<es_proc> ; es_name=<es_name> ; — 
N_fcycl=l ; NJcycl=1000 ; N_cylinc=l 
NSUB=0 ; NSTRT=0 ; NSJcycl=0 ; - 
NPRT=1) 

♦pack 1 
fxqt exit 
\endinput 


C2.1 Procedure to perform the monotonic loading calculations (file name pffdm.clp) 

♦procedure PFFDM ( es_proc ; es_name ; — 

N_fcycl ; N_lcycl ; N^cylinc ; — 

NSUB ; NSTRT ; NS Jcycl ; NPRT ) 

File to control monotonic loading to failure 

N_fcycl: First cycle number 

N_lcycl: last cycle number 

N_cylinc: cycle increment 

NSUB: subincrement flag (=0, to bypass) 

NSTRT: cycle to start subincrements(=0, to bypass) 

NS_lcycl: number of subincrements(=l, to bypass) 

NPRT: output storage cycle increment 

begin loop here 

♦set echo=off 

♦def icount = 0 

♦DO :CYCLOOP $NCYL = <[N_fcycl]>, <[N_lcycl]>, <[N_cylinc]> 

♦def icount = ( <icount> + 1 ) 

♦if < <icount> /eq <[NPRT]> > /then 
♦def iprint = 1 
♦def icount = 0 
♦else 

♦def iprint = 0 
♦endif 

*def/i SITCYCL = 1 
*def/i DUMIT = 1 
♦def/i $SNCYL = 1 


2 4 


TT 



*def delinc = <fN_cylinc]> 

♦def/d sLdisp == STAT.DISP* 

♦def/d ap.forc = APPL.FORC.* 

*def/a nst.disp = NST.DISP 
*def/a ost_disp == OST.DISP 
*def/a dst.disp == DST.DISP 
♦DO $ITCYCL= 1,1000 
♦def/i drfit = 1 
[xqt DRF.MNR 
select /PRINT = 0 
select /DSTATUS = 22222 
select /XFACTOR= 1.0 
select /DUM1T = <$ITCYCL> 
stop 

******************************************** 

. postpone iterations until damage occurs 

*G2M /name=drfit /type=D /maxn=l~ 

1, DRFITU.EX47.*, DRFDAT.l 
♦if < <drflt> /eq 0.0 > /then 
* def/i $ITCYCL = 1000 
♦endif 

******************************************** 

******************************************* 

iterate until change in damage resultant 
. forces are small 

*G2M /name=drfitn /type=I /maxn=l — 

1, DRFITV.EX47.*, DRFDAT.l 
♦print 1 DRFTTV.EX47.* DRFDAT.l /m=l 
♦if < <$ITCYCL> /ge 2 > /then 
*if < <drfitn> /eq 0 > /then 
♦def/i SITCYCL = 1000 
♦endif 
*endif 

******************************************** 

♦if < <$ITCYCL> /ge 2 > /then 
♦G2M /name=ddisp /type=D /maxn=l -- 
1, DST.DISP, DATA.l /m=49 
♦endif 

♦if < <$NCYL> /ge 2 > /then 
[xqt vec 

COMBINE <ost_disp> <- <st_disp> 

♦if < <$NCYL> /eq 2 > /then 
COMBINE <dst_disp> <- <st_disp> 

♦endif 

♦endif 

[xqt SSOL 

♦if < <drfit> /ne 0.0 > /then 
♦if < <$ITCYCL> /ge 2 > /then 
[xqt vec 

COMBINE <dst_disp> <- <st_disp> 

COMBINE <nst_disp> <- <dst_disp> + <ost_disp> 

COMBINE <st_disp> <- <nst_disp> 

♦endif 

♦endif 

♦call STRESS (direction=l; location= INTEG_PTS; print=<false> ) 
[xqt DGI 

select /PRINT = 0 

select /INC_SIZE = <delinc> 

select /N_C Y CLE = <$NCYL> 
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select /NINCR = <$SNCYL> 
select /DSTATUS = 22222 
stop 


*if < <$ITCYCL> ft q 1000 > /then 
♦if < <IPRINT> ft q 1 > /then 
♦if < <$NCYL> /ge 100 > /then 

♦copy 1, PLYDTM.<ES_NAME>.<$NCYL>.<$ITCYCL> = 1, PDAT.<ES_NAME>.* 

♦copy 1, DISPM.<ES_NAME>.<$NCYL>.<$ITCYCL> = 1, STAT.DISP.* 

♦endif 

♦endif 

♦endif 

♦ENDDO 

:CYCLOOP 

. *set echo=off 
♦end 

C3. Residual Strength Analysis Output 

Centroidal values of the ply stresses, strains, and internal state variables for each element can be found in 
the PDAT.xxxx data set (renamed as PLYDTM.EX47.xxxx.yyyy). Nodal displacements are located in the 
STAT.DISP.xxxx data set (renamed as DISPM.EX47.xxxx.yyyy). These data sets are located in the *.101 
output file and can be retrieved as separate individual files. An example of data output for one damaged 
element from one of the PDAT.xxxx data sets is given below. 


Record DAT_PLY,97 of dataset PLYDTM.EX47. 190.50 

1: 1 .9230D+05 6.8000D+02 7.2000D+03 9.8001D-03 1.0500D-02 3.1344D-02 
7: O.OOOOD+OO 1.3139D-02 2.1344D-02 2.0601D+05 6.8000D+02 7.2000D+03 
13: 1.0500D-02 9.8001D-03 -3.1344D-02 O.OOOOD+OO 1.2664D-02 -4.1344D-02 
19: 2.9008P+Q4 -Q656D+03 5.0370P+Q2 2.5822P-Q2 -5.521 8P-Q3 6.9959P-Q4 
25: 2.4230D-02 0.0O00D+00 0.Q000D+0Q -1.0801D+05 6.8000D+02 -5.0370D+02 
31: -5.5218D-03 2.5822D-02 -6.9958D-04 O.OOOOD+OO 2.3558D-02 O.OOOOD+OO 
37: 1.4610D+05 5.0053D+03 7.2000D+03 1.7683D-02 2.6170D-03 2.7494D-02 
43: 1.031 ID-02 1.3219D-03 1.7494D-02 1.5110D+05 1.8871D+03 7.2000D+03 
49: L8289D-02 2.01 1 ID-03 -2.6794D-02 1.0610D-02 3.0906D-03 -3.6794D-02 
55: 2.9008D+04 -6.8656D+03 5.0370D+02 2.5822D-02 -5.521 8D-03 6.9959D-04 
61: 2.4230D-02 O.OOOOD+OO O.OOOOD+OO 1.5110D+05 1.8871D+03 7.2000D+03 
67: 1.8289D-02 2.01 1 ID-03 -2.6794D-02 1.0610D-02 3.0906D-03 -3.6794D-02 
73: 1.4610D+05 5.0053D+03 7.2000D+03 1.7683D-02 2.6170D-03 2.7494D-02 
79: 1,031 ID-02 1.3219D-03 1.7494D-02 -1.0801 D+05 6.8000D+02 -5.0370D+02 
85: -5.5218D-03 2.5822D-02 -6.9958D04 O.OOOOD+OO 2.3558D-02 O.OOOOD+OO 
91: 2.9008D+04 -6.8656D+03 5.0370D+02 2.5822D-02 -5,521 8D-03 6.9959D-04 
97: 2.4230D-02 O.OOOOD+OO O.OOOOD+OO 2.0601D+05 6.8000D+02 7.2000D+03 
103: 1.0500D-02 9.8001D-03 -3.1344D-02 O.OOOOD+OO 1.2664D-02 -4.1344D-02 
109: 1.9230D+05 6.8000D+02 7.2000D+03 9.8001D-03 1.0500D-02 3.1344D-02 
115: O.OOOOD+OO 1.3139D-02 2.1344D-02 

This data is for element number 97 at load step number 190 and at the 50th (or last) iteration. Based on 
the structure of data sets in Appendix B Section 8, the underlined data is the output for the third ply (in this 
case that would be the 0° ply). The number 2.9008D+04 means that o n , the normal stress in the fiber 
direction, is 29,008 psi. This is 10% of the failure strength because we know that fiber fracture has occurred 
from the value of a n , 2.4230D-02. 

An example of the nodal displacements for the first nine nodes are listed below. Each row of data 
belongs to one node. The order of data is such that the x-direction translation is listed first, then the y and z 
translations followed by the rotations about the x, y, and z axes, respectively. 


Record DATA.l of dataset DISPM.EX47. 190.50 

1: 1.9339E-02 3.7506E-03 O.OOOOE+OO 0.0000E+00 0.0000E+00 0.0000E+00 


26 



7 

1.9377E-02 

3.2392E-03 

O.OOOOE+OO 

O.OOOOE+OO 

O.OOOOE+OO 

O.OOOOE+OO 

13 

1.9422E-02 

2.7315E-03 

0.0000E+00 

O.OOOOE+OO 

O.OOOOE+OO 

O.OOOOE+OO 

19 

1.9666E-02 

O.OOOOE+OO 

0.0000E+00 

O.OOOOE+OO 

O.OOOOE+OO 

O.OOOOE+OO 

25 

1.9475E-02 

2.2355E-03 

0.0000E+00 

O.OOOOE+OO 

O.OOOOE+OO 

O.OOOOE+OO 

31 

1.9656E-02 

4.2360E-04 

0.0000E+00 

O.OOOOE+OO 

O.OOOOE+OO 

O.OOOOE+OO 

37 

1.9532E-02 

1.7568E-03 

O.OOOOE+OO 

O.OOOOE+OO 

O.OOOOE+OO 

O.OOOOE+OO 

43 

1.9628E-02 

8.5391E-04 

O.OOOOE+OO 

O.OOOOE+OO 

O.OOOOE+OO 

O.OOOOE+OO 

49 

1.9585E-02 

1.2969E-03 

O.OOOOE+OO 

O.OOOOE+OO 

O.OOOOE+OO 

O.OOOOE+OO 


How often such data is stored in data sets is up to the user and is controlled by the NPRT variable in the 
runstream and the *copy 1 command in the procedure pffdm.clp . 
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Table 1. Material Properties of Unidirectional Ply of IM7/5260 


Ej^Msi 

E 22 , Msi 

G 12 , Msi 

V 12 

tpty. * n 

Pllcrit 

P22crit 

Growth law parameters: 
E 


22.16 

1.26 

0.75 

0.333 

0.006 

0.015 

0.008 

1.1695 


n .... 
dpara 


5.5109 
3.8686 x 10 7 


Table 2. Maximum Fatigue Loads Employed in Sample Calculations 


Layup 

Specimen 

geometry 

Maximum fatigue 
load (R = 0.1) 

[0/± 45/90] s 

unnotched 

3300 lb/in 

open hole 

2000 lb/in 

[0/90 3 ] s 

unnotched 

2480 lb/in 

open hole 

1572 lb/in 


Table 3. Lamina Material Properties 


Material 

En 

(Msi) 

E 22 

(Msi) 

G 12 

(Msi) 

v f 

(%) 

V 12 

90 
E c r 
(%) 

0 

Ecr 

(%) 

*° 

Yl2 

(%) 

* CT 
*Yl2 

(%) 

a AS4/8553-40 

19.7 

1.31 

0.65 

58.2 

0.34 

0.87 

1.56 

1.00 

10.00 

a AS4/938 

19.6 

1.36 

0.72 

57.2 

0.32 

0.50 

1.48 

1.00 

10.00 

b AS4/3501-6 

20.0 

1.36 

0.87 

60.3 

0.28 

0.50 

1.50 

1.00 

10.00 


a Eji, E 22 , G 12 , V 12 , e CT , and e cr measured by Boeing (ref. 16) 

b Eh, E 22 , G 12 , Vi 2 , £ct° » an< i £cr measured by Lagace et al. (ref. 17) 

o cr 

* Yl 2 and y 12 are approximations from Iosipescu shear test data (ref. 18) 





Input Finite Element Mesh, Current Damage State, 
Loading Condition, and Material Properties pj- 



Iterate Until 
damage resultant 
forces <0.1 


Calculate Ply-Level Strains, Stresses, 
and Damage Evolution 


Update Damage State 


No X LaminatX N o A SquilibriunK 
\ Failure ? / ' \ Established L 


Figure 1 . Progressive failure analysis scheme. 










(a) Dimensions and boundary conditions. 



X 


(b) The finite element model in sample calculation. 


Figure 2. Conditions and model of cross-ply laminated composite plate. All linear dimensions are in inches. 


32 


TIT 












y 

X 

Figure 6. Finite element model for a laminate with a central circular hole. 
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Number of cycles 


Figure 7. Stiffness loss of IM7/5260 laminates with central circular notch. 



Figure 8. Finite Element Mesh of the Center-Crack Tension Panel. 
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□ experimental 

□ model prediction 
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Figure 9. Residual Strength Predictions for the Center-Crack Tension Panels. 
10 6 

Stress, psi c 
10 5 

10 4 

10 3 
102 

0.001 0.01 0.10 1.00 10 100 

Strain, % 

Figure 10. Load Carrying Capability After Ply Failure. 
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